clc; clear; close all;
% ——— 参数设置 ———
x0      = 0.463442265;
y0      = 0.04532285;
z0      = 0.002136285;
b       = 30;
Nr      = 800;
r_list  = linspace(3.1,3.9,Nr);
Ntrans  = 500;
Niter   = 1000;

LEs = zeros(Nr,3);

parfor ir = 1:Nr
    r = r_list(ir);
    LEs(ir,:) = LogisticLEs(...
                    r, x0, y0, z0, b, Ntrans, Niter);
end

% ——— 绘图 ———
figure; hold on;
plot(r_list, LEs(:,1), 'r-', 'LineWidth',1);
plot(r_list, LEs(:,2), 'g-', 'LineWidth',1);
plot(r_list, LEs(:,3), 'b-', 'LineWidth',1);
xlabel('\gamma','FontSize',12);
ylabel('LEs','FontSize',12);
legend('LE_1','LE_2','LE_3','Location','Best');
box on; axis tight;


function LE = LogisticLEs(r, x0, y0, z0, b, Ntrans, Niter)

    % ——— 丢弃暂态 ———
    x = x0; y = y0; z = z0;
    x_prev = x; z_prev = z;
    for i = 1:Ntrans
        [x_new,y_new,z_new] = Logistic(x,y,z,x_prev,z_prev,r,b);
        x_prev = x; z_prev = z;
        x = x_new; y = y_new; z = z_new;
    end

    % ——— QR 方法初始化 ———
    Q      = eye(3);
    sumLog = zeros(1,3);

    % ——— 主循环：每步迭代 + QR ———
    for i = 1:Niter
        % 1) 系统迭代
        [x_new,y_new,z_new] = Logistic(x,y,z,x_prev,z_prev,r,b);

        % 2) 直接用解析雅可比
        J = LogisticJaco(x, y, z, x_prev, z_prev, r, b);

        % 3) 切空间传播并正交化
        [Q,R] = qr(J * Q, 0);
        sumLog = sumLog + log(abs(diag(R)))';

        % 4) 更新状态
        x_prev = x;   z_prev = z;
        x = x_new;    y = y_new;    z = z_new;
    end

    % ——— 平均化得到 LE ———
    LE = sumLog / Niter;
end

function J = LogisticJaco(x,y,z,x_prev,z_prev,r,b)
    c1 = exp(-2*b);
    c2 = exp(-b)*r;
    xdot = x - x_prev;
    J11 = r*(1-2*x);    J12 = -r;        J13 = 0;
    J21 = c2*(-2*y -2*z + z_prev);
    J22 = -c1 + c2*(2 -2*x + x_prev);
    J23 = c2*(-2*x + x_prev);
    J31 = c2*(-2*z^2 -2*y -1);
    J32 = c2*(-2*x);
    J33 = -c1 + c2*(2 -4*x*z +4*x_prev*z);
    J = [J11, J12, J13; J21, J22, J23; J31, J32, J33];
end
